Homogeneous nucleation in supersaturated vapors of methane, ethane, and 
carbon dioxide predicted by brute force molecular dynamics 

Martin Horsch'^, Jadran Vrabec*^, Martin Bernreuther-'-, Sebastian Grottel^, 
Guido Reina^, Andrea Wix', Karlheinz Schaber', and Hans Hasse^ 



o 
o 

<N 

Oh' 
< 



^ ■ 



Abstract 

Molecular dynamics (MD) simulation is applied to the condensation process of supersaturated vapors of methane, 
ethane, and carbon dioxide. Simulations of systems with up to a million particles were conducted with a massively 
parallel MD program. This leads to reliable statistics and makes nucleation rates down to the order of 10^" m~^s~^ 
accessible to the direct simulation approach. Simulation results are compared to the classical nucleation theory (CNT) 
as well as the theory of Laaksonen, Ford, and Kulmala (LFK) which introduces a size dependence of the specific surface 
energy. CNT describes the nucleation of ethane and carbon dioxide excellently over the entire studied temperature 
range, whereas LFK provides a better approach to methane at low temperatures. 
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Homogeneous nucleation was discussed theoretically by 
Gibbs [Ij and studied in depth by Volmer and Weber [2j as 
well as Farkas [3 . In combination with experiments car- 
ried out by Wilson [4J and Powell [5j during the same pe- 
riod, these efforts established the classical nucleation the- 
ory (CNT), which is known to be accurate in many cases 
but to fail in others (H [TJ [8] . 

Molecular simulations are applied to this problem since 
the late 1950s, when Alder and Wainwright [9j observed 
a first-order phase transition in molecular dynamics (MD) 
simulations of the hard sphere fluid. In the 1970s, McGinty 
[lOj studied liquid clusters of the Lennard- Jones (LJ) fluid 
in MD simulations, and Rao et al. [11] described the con- 
densation of a supersaturated vapor with results obtained 
from both Monte Carlo (MC) and MD simulations. Some 
common approaches to the dynamics of nucleation, such as 
MD simulations with inserted droplets [T2l [T3l [Til [T5l [TB] 
'pi transition path sampling [17^ 18^ 19] as well as MC sim- 
■ulations [H [20l [HI [22 [23] , do not lead immediately to 
'the velocity of the phase transition, but only to indirect 
information, e.g. on the required activation energy. The 
present study discusses brute force MD simulations, which 
are aimed at the direct reproduction of a nucleation pro- 
cess by means of the deterministic simulation of a large 
system. 

Nucleation processes are characterized by the nucle- 
ation rate J, i.e. the number of stable liquid nuclei gener- 
ated per volume and time, and their critical size i.e. the 
number of molecules in a nucleus with maximal Gibbs en- 
ergy of formation. Droplets above that size have a higher 
probability to further grow, whereas smaller clusters tend 
to evaporate. Due to current limitations in the available 



computational resources, only nucleation processes with 
extremely large values of J can be simulated by MD. How- 
ever, when nucleation occurs very rapidly, the vapor phase 
is not fully in equilibrium with the emerging droplets and 
the critical size is not constant. It is nonetheless possi- 
ble to determine nucleation rates if one follows the some- 
what heuristic approach proposed by Yasuoka and Mat- 
sumoto |24j. Most recent direct MD studies of nucleation 
[23 ESI EH [^ adhere to this method. 

The method of Yasuoka and Matsumoto requires sys- 
tem sizes and simulation times to be as large as possible. 
Due to restrictions of computational power, the lowest nu- 
cleation rates which can be obtained nowadays with this 
approach - above 10^^/(m^s) in the present study - ex- 
ceed those which actually can be observed in experiments 
by about seven orders of magnitude [29] . This gap can only 
be closed by predictions on the basis of nucleation theories 
that express the dependence of J and ^* on temperature 
and pressure, where the latter is often given in terms of 
the supersaturation Sp{p^ T) = p/pcr{T)^ i.e. related to the 
vapor pressure Pcr- Reviews following the progress of the 
last decades are provided by Oxtoby [30l [31] and Ford [8] . 
For a description of advanced experimental methods see 
also Fladerer and Strey [32j as well as Hand [29j. 

In the following sections, CNT and a version of the 
Dillmann-Meier [33j model due to Laaksonen, Ford, and 
Kulmala [34j, referred to as LFK here, will be discussed 
and compared to data from direct MD simulations. It is 
also necessary to comment on the mean first passage time 
(MFPT) approach, an indirect method which fits a pre- 
defined kinetic model with three parameters to simulation 
results [35j [36] . 
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Figure 1: left - Gibbs energy of cluster formation in CO2 at 253 K for pressure values of 2100, 2350, and 3350 kPa 
according to CNT ( — ) and LFK ( — ); center - Strong (256000 molecules) and weak scaling (2000 molecules per 
process) of Isi simulating CH4 at 102 K and 0.730 mol/1 on the Intel Xeon cluster mozart at the chairs Simulation of Large 
Systems and Numerics for Supercomputers^ Universitat Stuttgart; right - Strong (256000 molecules) and weak scaling 
(4000 molecules per process) of Isi for 2CLJQ fluid models on the Intel Xeon cluster cacau at the High Performance 
Computing Center Stuttgart: C2H6 at 183 K and 0.365 mol/1 ( — ) as well as CO2 at 253 K and 1.670 mol/1 ( • • • ) 



2 Nucleation theories 



2.1 Classical nucleation theory 

The foundations of CNT were laid by Gibbs [Ij and fur- 
ther developed by Volmer and Weber [2]. Important sub- 
sequent contributions were made by Farkas [3] , Becker and 
Doring [37J, Zel'dovich [38j, and Feder et al. For the 
further development of the theory compare Kashchiev [40] 
and Vehkamaki [41j. 

The starting point of this theory is the capillarity ap- 
proximation: the dispersed liquid phase, composed of the 
clusters emerging during nucleation, is assumed to have 
the same thermodynamic properties as the saturated bulk 
liquid. It is also assumed that all liquid clusters are spher- 
ical. CNT describes how, under such preconditions, nu- 
cleation rate and critical size depend on temperature, su- 
persaturation, and a few properties of the fluid which are 
independent of 5p, such as the planar interface surface 
tension 70 and the densities pa- and pi, referring to the 
saturated vapor and liquid, repsectively. 

The surface energy <P{l) of a cluster with l molecules, 
the surface area A{l)^ and the specific surface energy £ 
amounts to £A{l). The capillarity approximation assigns 
f = 70, and CNT further assumes spherical clusters, hence 
A{l) = ^^/tt (6^/pi)^/^. The Gibbs energy of cluster forma- 
tion in a supersaturated vapor consists of a negative bulk 
contribution and a positive surface contribution [Tj ^ . It 
amounts to 



AG, 



W-^(i) + (i-OAm, 



(1) 



and reaches a maximum at the size of the critical nu- 
cleus. It can be seen from Figure [T] (left) that the criti- 
cal size is strongly dependent on the supersaturated vapor 
pressure; it diverges as the supersaturated vapor pressure 
p approaches the saturated vapor pressure Pa- of the bulk 
fiuid. The chemical potential difference. 



/ dp/ P., 



(2) 



is an integral between pa- and p at constant temperature. 
In metastable equilibrium, the ^-cluster number density 



= NJV^ where is the number of clusters with ex- 
actly i molecules, amounts to 



pi exp 



-AG, 
kT 



where pi can be estimated from 



E 



(3) 



(4) 



The impact rate 13 of vapor molecules on a cluster per 
surface area can be approximated by 



p 



V2TrmkT' 



(5) 



where m is the molecular mass [42 . Assuming that every 
collision of a monomer with a critical nucleus leads to the 
formation of a cluster with 1 + ^* molecules, the nucleation 
rate is given by 

J = p*pA*Z^. (6) 

Here and elsewhere, all quantities marked with an aster- 
isk refer to critical nuclei. The factor [5 A* expresses the 
impact frequency of monomers on a critical nucleus, or 
equivalently, the rate at which critical nuclei grow to a 
size of 1 + ^* molecules. 

The two remaining factors, Z and "d^ represent correc- 
tions with respect to the nucleus density, the kinetics of 
the nucleation process, and the temperature inside a nu- 
cleus. The metastable equilibrium breaks down near the 
critical size, and the actual number density of critical nu- 
clei is considerably lower than their metastable equilibrium 
density p* . The Zel'dovich factor. 
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takes into account both this difference and the probability 
that that a nucleus above the critical size does not continue 
to grow [38j. 
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Nuclei that reach the critical size usually have grown 
very fast and retain part of the latent heat. Let Cy be the 
isochoric heat capacity of the vapor and AHy the enthalpy 
of evaporation. From considerations of Feder et al. [39j it 
follows that on average, critical nuclei are overheated by 



AH, ' 



(8) 



This increase in temperature reduces the nucleation rate, 
an effect which is quantified by the non-isothermal factor 
The energy added to a critical nucleus when it acquires 
an additional molecule is 



q = AH^ - — 



d<p{i) 
dL 



in excess of what is needed to extend its area at the same 
temperature. The standard deviation of the energy of va- 
por molecules that collide with a cluster is 



Finally, the non-isothermal factor is given by 



62 + g2 • 



(10) 



(11) 



2.2 Model proposed by Laaksonen, Ford, 
and Kulmala 

The LFK model ^34j is a version of the Dillmann-Meier ap- 
proach [33] which postulates the surface energy of a cluster 
with ^ G N molecules to be 



{i) = Jo A{i) -\- rkT Int. 



(12) 



The adjustable parameters of this model are r and hz{L) for 
^ G N, as well as pi which is expressed indirectly by means 
of a normalization parameter qq. By comparing the Fisher 
j43] equation of state. 



P 



(13) 



to a virial-type expansion of second order values for i<i{l) 
and i<i{2) are defined. Laaksonen et al. [34J represent this 
in terms of the monomer density as 



Pi 



kT 
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where the second virial coefficient is, in this case, defined 
as 



B 



p-^kT. 



(15) 



They obtain the expressions 
'^(2) = (=^-'^(1)0 
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with = 7oA(l)/fcT, by applying an approximation to 
Equation This is extended to higher k{l) by 



K{i) 



1 + aii 



-1/3 , 



-2/3 



(16) 



which are the first three contributions of an expansion in 
terms of the inverse cluster radius ~ The coef- 

ficients ai and a2 are determined by equating the expres- 
sions for and hz{2) with Equation (p!6|) . 

Note that since Equation ([12]) multiplies with 
A{l) ~ the value of 0^2 only infiuences a constant 

summand which cancels out in the expression for AG^. 
Laaksonen et al. [34] proposed r = 0, and Ford et al. 
[44] showed that with this particular assignment, the pa- 
rameter go cancels out as well. To obtain a convenient 
expression, we set qo = p^/kT^ which leads to 



^(1) 

/t(2) 



70^(1)' 

-2Bp^ - kT\n{-Bp^/kT) 



The Zel'dovich factor takes the form 



Z = — 

3^* 
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(17) 

(18) 
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and the energy released on addition of a monomer to a 
cluster amounts to 



AH„ - — 



kT 
2" 



70^(1) 



(20) 



Compared to CNT, the LFK model hence introduces a size 
dependence of the specific surface energy which is governed 
by the single parameter ai. Figure [T]( left) illustrates that 
this size dependence becomes particularly significant at 
high supersaturations, where the critical nucleus is small. 

2.3 Mean first passage times 

Let us next consider the kinetics of a nucleation process. 
For a supersaturated fluid in a volume V that exhibits the 
nucleation rate J, it might be expected that the first stable 
nuclei appear on average after a temporal delay, expressed 

by 

r(0 « jy, (21) 

for some l > ^ after the onset of nucleation [45j. The av- 
erage delay F {i) until the first cluster with t molecules ap- 
pears is called the mean first passage time of i. Wedekind 
et al. [36^ 46, 47J generalized this approach to a theory of 
condensation processes, here referred to as MFPT. Bartell 
and Wu [35j obtained an identical result for freezing, and 
Zhang et al. [48j applied it to melting processes. According 
to this approach, the mean first passage time is approxi- 
mated using a Gauss error function. 



r(0 



;i + erf(x(^-^*))]. 



In particular, this approach leads to 



and thus 



lim r{i) 



J 



2F" 



1 



(22) 



(23) 



(24) 



2r*F' 

from Equation ([2T]) with ^ ^ 00 [45j. These approxima- 
tions are intended to hold only 'in the vicinity of the crit- 
ical size' and '[ujnder reasonably high barriers' [36j. 



3 





model 


m [u] 


r A 1 

a [A] 


e [A: X K] 


Q [B] 


7- TAT 

L [A] 


CH4 


LJ 


16.04 


3.7281 


148.55 






C2H6 


2CLJQ 


2 X 15.03 


3.4896 


136.99 


0.8277 


2.3762 


vyW2 


9PT TO 


9 V 99 nn 




1 '^'^ 99 
±00. zz 


Q 7QQQ 

0. / yoo 


9 Ai 



Table 1: Parameters of the molecular models for methane, ethane, and carbon dioxide 




Figure 2: left - Number of nuclei containing at least 75 molecules in supersaturated CO2 vapor over simulation time; 
center - Number of nuclei containing at least 25, 100, . . . , 1000 molecules over simulation time in (63.7 nm)^ filled with 
methane at 130 K and 1.606 mol/1; right - Cluster formation delay i^xit) for CO2 at 253 K and 3.150 mol/1 with 
AT = 884700 and X G {2-^,2-^,2-^^}. 



3 Simulation method 



Methane, ethane and carbon dioxide were selected in 
the present work because of their qualitatively different 
molecular properties. Methane is almost spherical and 
weakly octupolar, thus it can be described by a single 
Lennard- Jones site with the pair potential 

Ethane molecules are dumbbell-shaped and thus sig- 
nificantly anisotropic in geometry but only slightly 
quadrupolar. Carbon dioxide molecules are both strongly 
anisotropic and quadrupolar. The inter molecular interac- 
tions of these two fiuids were described by the two-center 
LJ model with an embedded point quadrupole (2CLJQ). 
Additional parameters of the 2CLJQ model are the molec- 
ular elongation L and the quadrupole moment Q. The 
parameters of the molecular models, cf. Table [TJ were ad- 
justed to experimental vapor-liquid equilibria in prior work 

Series of MD simulations of nucleating vapors were con- 
ducted using a version of the Isi program [50|. The sim- 
ulations were carried out in the canonical ensemble, with 
a time step between 3 and 7 fs, depending on the system 
temperature. The cutoff radius was larger than 4.5cr in 
all simulations. The temperature of the whole system was 
kept constant by isokinetic scaling. 

To follow the kinetics of the phase transition in de- 
tail, a criterion which detects clusters of molecules, i.e. 
the dispersed liquid phase, must be applied to the whole 
ensemble. In past studies, a considerable number of 
different cluster criteria were discussed and compared 
I5ll[52l[53l[5l[^. Those presented by Hih |56j and Still- 



inger [57] are among the most common ones. They are 
applied to all pair interactions and the clusters are de- 
fined as the connected components of the graph with the 
molecules as its nodes and edges between the pairs with 
interactions that fulfill a pair critierion. The Hill energetic 
criterion is defined by 

o 

TYlll 

<ri,) + ^ < 0, (26) 
and the Stillinger geometric criterion by 

nj < rgc, (27) 

where rgc = 1.5<j for the Lennard- Jones fiuid. The above 
definitions distinguish the bulk phases. A hybrid cluster 
criterion which combines these definitions was consistently 
observed to select only few clusters with extremely short 
lifetimes, whereas it reliably detected stable clusters of all 
sizes. It is defined as follows: 

• All molecules i for which the energetic single- 
molecule criterion 

mvl ^T.i^jU^,,{rij) < 0, (28) 

holds are defined to be liquid. 

• Two liquid molecules are regarded as connected 
whenever they fulfill the Stillinger criterion. For 
the 2CLJQ model the maximal connection radius is 
given by 
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• Clusters are determined by covering a graph con- 
sisting of these connections with maximal bicon- 
nected components and eliminating their overlap. 
Monomeric clusters are regarded as vapor molecules. 

Biconnected components are, by definition, subsets of a 
graph that cannot be separated into two unconnected parts 
by removing only one vertex. This reflects the idea that a 
cluster should not consist of several sub-nuclei connected 
by a single molecule, because structures that do depend 
on such a connection tend to be extremely unstable. 

The Hill energetic criterion favors molecules with a low 
kinetic energy, and hence leads to artefacts in the cluster 
temperature, i.e. clusters are observed to be colder than 
they actually are. This effect carries over to the hybrid 
criterion. For this reason, temperature data as displayed 



in Figure [3] (center) were gathered by applying only the 
geometric and the biconnectivity parts of the hybrid cri- 
terion. 

The MD program Isi relies on spatial domain decom- 
position for parallel simulations [50] . The operation of par- 
titioning a very large graph into biconnected components 
was handled by including the boost library [58j which im- 
plements Tarjan's linear time algorithm. Figure [1] (center 
and right) shows that Isi, both with and without cluster 
recognition, scales well on typical clusters of workstations. 
This permits simulations of volumes V ~ 10~^^ m^ for a 
time span of a few nanoseconds with an acceptable com- 
putational effort. Thus with the direct approach, which 
requires at least some stable nuclei to appear, only val- 
ues J > 10^^/(m^s) are accessible unless correspondingly 
larger computational resources are employed. 




150 200 250 300 50 100 150 200 250 3100 3200 3300 3400 
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Figure 3: left - Cluster net growth rate as a function of cluster size and time in C2H6 at 280 K and 2.800 mol/1; center 
- Cluster temperature as a function of cluster size and time in C2H6 at 280 K and 2.800 mol/1; data (circles) and running 
average ( — ) from 83 to 250 ps as well as data (squares) and running average ( — ) from 650 to 700 ps after the initial 
state; right - Nucleation rate of C2H6 at 280 K; small squares: J(50), large squares: J(75) and J(IOO) 



4 Simulation results 



4.1 Growth rates of single nuclei 

Both CNT and the Dillmann-Meier model assume that the 
properties of clusters with a given size depend only on the 
temperature and supersaturation of the vapor. Hence, one 
should expect droplets of the same size generated earlier 
and later during a simulation run to have, on average, con- 
stant temperatures and rates of growth and evaporation. 

It is known from MD simulations by Tanumura et al. 
[59] that this does not necessarily hold in the initial period: 
the very first clusters of a given size have a higher kinetic 
energy than those which belong to the actual metastable 
vapor. This is due to the fact that molecules lose potential 
energy when they attach to a droplet, which transforms to 
kinetic energy. The first large clusters have experienced a 
relatively fast growth process and hence retain more of this 
latent heat. The present simulations confirm this observa- 
tion, cf. Figure [3] (center), which shows that the largest 
existing clusters have a temperature of 282 K and above, 
while the whole system temperature is fixed at 280 K. The 
lower curve, collected between 650 and 700 after the simu- 
lation onset, exhibits a local maximum of the cluster tem- 
perature at a size of about 70 molecules. The temperature 
of smaller clusters changes comparatively little over simu- 



lation time, i.e. with respect to the higher curve, whereas 
for larger clusters it decreases considerably. The tempera- 
ture of the clusters with ^ < 70 has reached a steady state, 
but no thermal equilibrium with the vapor, which implies 
that these clusters are unstable. For l > 70, no steady 
state is established and the cluster temperature approaches 
280 K. Hence, the local maximum of the temperature plot 
marks the transition between unstable clusters and stable 
nano-droplets, i.e. the size of the critical nucleus. It agrees 
well with the critical size of 78 indicated by CNT, as op- 
posed to the LFK model which predicts = 207, cf. Table 

m 

As an effect which is closely related to the overheating 
of the dispersed phase, clusters generated early in the con- 
densation process evaporate at a higher rate than those 
which are generated later, cf. Figure [3] (left). These data 
were collected during the same simulation as those from 
Figure [3] (center). Note how the growth rate of clusters 
is negative for sizes significantly larger than the critical 
nucleus, where l* can be estimated either from the tem- 
perature profile any of the theories. This phenomenon 
was also observed by Yasuoka and Matsumoto [24]. The 
positive contribution to cluster growth, which is due to 
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Table 2: Comparison of the nuclation rate from an MFPT indirect analysis according to Wedekind et al. [46j with the 
value J(10) from another simulation of that system analyzed with the method of Yasuoka and Matsumoto; a single set of 
values for the LJ fluid is interpreted as both argon and methane 
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Table 3: Simulation results and theoretical values of nucleation rates for supersaturated methane at 170 K with pcr= 2328 
kPa, pcr= 2.429 mol/1, and 70= 2.07 g/s^; bold values: threshold l > according to theory 



condensation of the vapor phase, remains constant over 
simulation time. This indicates that while the tempera- 
ture of clusters with a given size changes, the temporal 
evolution of the system does not significantly affect their 
surface area. From the decrease in cluster evaporation over 
time, it clearly follows that for MD simulations starting 
from a cluster-free configuration, the first clusters differ 
significantly from the much larger number of clusters of 
the same size that appear at a later stage of the process. 

Changing rates of evaporation also imply that the crit- 
ical size can actually vary during a simulation run with 
a very high supersaturation. For systems at very high 
supersaturations, it is impossible to observe a metastable 
vapor phase, because nucleation begins immediately. Such 
phenomena are only realistic if it is technically possible to 
increase the supersaturation faster than the vapor phase 
can produce small clusters. 

A magic number effect can be observed for small clus- 
ters of both the LJ and the 2CLJQ fiuids: the rate of 
evaporation is comparatively low for clusters with t G 
{8, 11, 14, ... , 26} molecules. As opposed to this, clusters 
with i G {4,9,12,15} are detected to have particularly 
high rates of evaporation. Within this range, 23 and 26, 
but also 13 and 19, are known as preferred cluster sizes of 
the L J fluid [6OJ . In the present study, the magic numbers 
may well be a side effect of the biconnectivity requirement 
of the hybrid cluster criterion. This conclusion is also sug- 
gested by the fact that the observed magic numbers do not 
depend on the employed molecular model. 

4.2 Nucleation rates 

A nucleation process at constant pressure in an infinitely 
large system occurs, by definition, with the nucleation rate 



J = lim lim 



dt 



(30) 



From molecular simulation in the canonical ensemble, a 
smoothed where the statistical noise is reduced, can 

be constructed from by averaging over a number of 
time steps. Such smoothed plots are shown in Figure [2] 



(left and center). The nucleation rate may then be ap- 
proximated by the expression 



j(0 



max — ; — . 

*>*o at 



(31) 



This approach was introduced by Yasuoka and Matsumoto 
[24j . The values of J(^) are meaningful for all i> i* . How- 
ever, it has to be taken into account that as the condensa- 
tion proceeds in a closed system, the number of monomers 
decreases and the pressure in the vapor is reduced signifi- 
cantly, which causes larger nuclei to be formed at a lower 
rate, cf. Figure [2] (center). 

The present simulation data suggest that, as expected, 
the values of J{i) are similar for values of i above a certain 
value, except for cases where the supersaturation decreases 
significantly, cf. Figures [3] (right) and [H (left and center) as 
well as Tables [31 m O and[6l On the other hand, J{i) with 
very small i is often significantly elevated, which raises 
doubts whether results related to J(6), cf. Kraska [27] . 
can lead to reliable conclusions. The nucleation rates are 
displayed together with pressure values, which were taken 
in the center of the interval where the plot of pi{t) and the 
linear approximation from which the value of J[i) is ob- 
tained, roughly agree - for instance, after two nanoseconds 
in the case of Figure [2] (left). 

Wedekind et al. [36j propose an indirect method for the 
determination of both the nucleation rate and the criti- 
cal size from simulation data on mean first passage times. 
This approach consists in fitting the values of Too, and 



such that Equation ([22]) agrees optimally with the actual 
plot of r{i) from an MD simulation. However, Wedekind 
et al. [47l Figure 4] also note that the size of the critical nu- 
cleus determined according to this MFPT based approach 
can deviate by a factor larger than two from the 'nucle- 
ation theorem,' 



d\nJ 
dhiSr) 



^*+2, 



(32) 



obtained by Oxtoby and Kashchiev [61 J in a similar ver- 
sion. That is not necessarily an argument against the 
method, since the nucleation theorem is known to be valid 
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2.352 


1.70 


1.47 


25 


6.7 X 10^2 


44 


6.9 X 10^2 


59 


2.5 X 10^1 


2.352 


1.69 


1.46 


100 


1.2 X 10^2 


46 


5.9 X 10^2 


61 


1.7 X 10^^ 


2.481 


1.72 


1.48 


25 


1.4 X 10^^ 


42 


9.1 X 10^2 


55 


6.1 X 10^^ 


2.481 


1.70 


1.47 


100 


2.9 X 10^2 


44 


6.9 X 10^2 


59 


2.5 X 10^^ 


2.609 


1.73 


1.48 


25 


2.2 X 10^^ 


41 


1.0 X 10^^ 


54 


6.7 X 10^^ 


2.609 


1.70 


1.47 


100 


7.6 X 10^2 


44 


6.9 X 10^2 


59 


2.5 X 10^^ 


2.695 


1.71 


1.47 


200 


8.4 X 10^2 


43 


8.0 X 10^2 


57 


3.8 X 10^^ 




1 70 


1 zL7 






AA 






Zi.O X ±u 



Table 4: Simulation results and theoretical values of nucleation rates for supersaturated carbon dioxide at 253 K with p^ = 
1961 kPa, p(j= 1.169 mol/1, and 70= 86.2 g/s^; bold values: threshold l > l* according to theory 



only for moderate supersaturations [62]. Table [3] com- 
pares a new MD simulation, evaluated according to the 
method of Yasuoka and Matsumoto, with data obtained 
by Wedekind et al. [46j following the MFPT approach. 
The value of J(10) is probably larger than J, since nucle- 
ation rates of about 10^^ - 10^^ usually imply critical sizes 
L* ^ 10. However, J(10) is significantly lower than the 
MFPT extrapolation. 

4.3 Delay of cluster formation 

Statistics on the formation delay of ^-clusters are shown in 
Figure [2] (right). The plots are of the type 

uS) = max{nGN| S(,>n)^p. > ^p}, (33) 

with < ^ < 1, i.e. they show the nucleation threshold 
z^^(t) passed by a mole fraction ^ of the condensing fluid 
at the time t. 

For instance, after 0.5 ns, A^/16 = 55300 or more 
molecules are in clusters with a size i > i^i/ie (0.5ns) = 450, 
but the clusters of 451 or more molecules contain less 
than 55300 molecules. At the same time, the thresh- 
old corresponding to A^/1024 = 864 molecules lies at 
z^i/1024 (0.5ns) = 1519, i.e. there are at least 864 molecules 



in clusters with l > 1519, but not in clusters with l > 1520. 
Thus the plot corresponding to i^i/io24(^) shows the devel- 
opment of the largest cluster. For that reason, it oscillates 
more than the other plots. 

For ^ ^ 0, the expected values of i^^{t) converge 
by definition towards the inverse function of F(^), since 
by inverting such a plot the first passage times are ob- 
tained. Consider such first passage times from simulations 
of methane, cf. Figure [4] (right). The data correspond to 
droplets much larger than the critical nucleus, cf. Tabled 
From Equation ([23|) and ([24|) it would be expected that 
the mean first passage time converges according to 

lim r(0 = (34) 

which corresponds to 66 ps for 4.116 mol/1 and to 38 ps for 
4.298 mol/1, if we accept the values of J(225) determined 
with the method of Yasuoka and Matsumoto. However, no 
convergence on such a timescale can be observed, and this 
is certainly not a matter of the statistical uncertainty of 
conducting a single simulation. The tendency of the mean 
first passage time to diverge instead of reaching a plateau 
can also be observed for data published by Wedekind [45l 
Figure 4.11 (bottom)]. 



5 Comparison of theory and simulation 



The simulation results for J{i) are compared to CNT 
and LFK in Figures [3] (right) and [H (left and center) as 
well as Tables O H [5l and El 

The correlation between p, p and A/i at constant T, 
which is necessary to evaluate the considered models, was 
obtained from simulations of small supersaturated systems 
analogous to those described in previous work [63]. The 
dependence of p on p between data points was approxi- 
mated by a linear fit. The resulting isotherms were used to 
estimate the density p(p, T) of the vapor, which decreases 
over simulation time, to reflect that with a decreasing su- 
persaturation, nuclei should be expected to emerge at a 
lower rate - note that the values of p shown in the ta- 
bles correspond to the density of the entire system, not to 
the remaining vapor. The isotherms were also applied to 
determine the second virial coefficient for the LFK model 



according to Equation ([151), chemical potential 

difference between the saturated and the supersaturated 
vapor according to Equation In Tables [SKU and [5l the 
supersaturation with respect to the density, Sp = p/ pa{T), 
as well as the pressure are shown together with 

f p{p,T)-p^{T) \ 
^/i = exp ( — 1 , (35) 

i.e. the supersaturation with respect to the chemical po- 
tential, where Pa{T) and Pa{T) refer to the saturated va- 
por at the given temperature. Occasionally, the identity 
Sp = Sp = Sp is assumed in the literature [H [HI |45] ; in 
particular, it is used for the derivation of Equation (|32]) . 
where Sp replaces the more accurate Sp. However, near 
the spinodal this is always a bad approximation, since dp/ 
dp ^ holds there by definition. 
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i 


T / \ r ^ 1 1 

J{i) [m '^s ^] 


^*(CNT) 


T / N TV T m \ r 'A. 1 1 

J(CNT) [m "=^8 ^] 


^*(LFK) 


J(LFK) [m ^=^8 ^] 


1.079 


1.077 


1.047 


25 


3.3 X 10^^ 


629 


2.3 X 10^^ 


1390 


2.7 X 10^^ 


1.079 


1.078 


1.047 


50 


6.5 X 10^2 


607 


3.1 X 10^^ 


1350 


5.7 X 10^^ 


1.177 


1.112 


1.067 


50 


4.8 X 10^2 


222 


1.5 X 10^2 


596 


2.4 X 10^^ 


1.177 


1.111 


1.066 


75 


1.1 X 10^2 


228 


1.4 X 10^2 


608 


1.7 X 10^^ 


1.238 


1.133 


1.078 


50 


2.4 X 10^^ 


140 


7.7 X 10^2 


416 


6.2 X 10^^ 


i .Zoo 


i . iOO 


1 nvQ 


/ 


1 A V 1 0^3 


loo 


o.O X iU 


^Uo 


1 n V 1 0^6 



Table 5: Simulation results and theoretical values of nucleation rates for supersaturated carbon dioxide at 285 K with p^- 
4712 kPa, pa= 3.270 mol/1, and 70= 24.5 g/s^; for all values the threshold l is lower than ^* according to theory 



T[K] p 



mol 



1 



p [kPa 



m3s 

10^ 

10^1 

10^1 

10^1 

lO^i 

1032 
1033 
1033 

10^ 
10^1 

1032 
1032 
1032 
1033 
1033 
1033 



(,*(CNT) J(CNT) 



10^ 

1030 
1030 
1030 
1030 

1032 
1032 

1030 
1030 

1032 
1032 
1032 
1033 
1033 
1033 



'(LFK) J(LFK) 



m3s 
10^ 

1Q32 

1Q32 
1Q32 
1Q32 

10^ 

1Q26 
1Q26 

10^7 

1Q28 
1029 
1029 



CH4 



106.0 
114.0 
114.0 
114.0 
114.0 
130.0 
130.0 
130.0 
130.0 



0.758 
0.851 
0.851 
0.925 
0.925 
1.432 
1.606 
1.693 
1.780 



503 
616 
614 
641 
629 
1022 
1095 
1148 
1167 



25 
75 
150 
75 
150 
700 
75 
25 
25 



1.8 X 

2.7 X 

2.8 X 
5.8 X 
5.5 X 

2.1 X 

6.2 X 
2.5 X 
4.0 X 



22 
23 
23 
22 
23 
31 
26 
24 
23 



1.6 X 

2.1 X 

2.0 X 
4.5 X 

3.2 X 

3.2 X 

1.3 X 
3.2 X 

4.1 X 



16 
19 
19 
18 
18 
31 
26 
22 
21 



9.9 X 
5.7 X 

5.4 X 
1.2 X 

8.5 X 

3.1 X 

1.7 X 

6.2 X 

7.8 X 



176.5 
176.5 
280.0 
280.0 
280.0 
280.0 
280.0 
280.0 



0.385 
0.400 
2.470 
2.470 
2.550 
2.800 
2.950 
2.950 



455 
467 
3283 
3278 
3307 
3397 
3430 
3427 



25 
25 
50 
75 
100 
100 
75 
100 



1.3 X 
2.0 X 
3.9 X 
1.7 X 

1.7 X 

1.8 X 
3.3 X 

1.9 X 



16 
16 
131 
135 
116 
78 
68 
69 



4.5 X 

6.4 X 

4.6 X 
4.2 X 

7.7 X 

2.5 X 

3.6 X 
3.5 X 



10^ 

10^1 

10^1 

1Q32 
1Q33 

1033 
1032 
1033 
1Q33 



10^ 

10^1 

10^1 

10^1 

1Q32 

1032 
1032 
1033 
1Q33 



15 
14 
322 
329 
290 
207 
186 
188 



2.7 X 
3.9 X 

2.8 X 
2.0 X 
1.2 X 

9.0 X 

3.1 X 
2.8 X 



10^ 

1Q30 
1Q30 
1Q30 

10^1 

1028 

1028 
1030 
1030 



CO2 



237.0 
237.0 
237.0 
237.0 
237.0 
269.0 
269.0 
269.0 
269.0 



1.700 
1.750 
1.850 
2.000 
2.450 
3.120 
3.120 
3.800 
3.800 



2283 
2317 
2322 
2333 
2499 
4142 
4131 
4350 
4343 



75 
75 
300 
300 
75 
25 
75 
50 
75 



1.1 X 
1.8 X 

9.6 X 
2.3 X 

4.7 X 
4.7 X 
3.6 X 
6.5 X 
4.1 X 



50 
48 
47 
46 
37 
81 
83 
54 
55 



9.4 X 

1.6 X 

1.7 X 
2.0 X 

1.4 X 
4.7 X 
4.2 X 
2.6 X 

2.5 X 



57 
55 
54 
53 
40 
147 
150 
99 
100 



7.7 X 
1.3 X 
1.5 X 

2.0 X 
4.2 X 
6.7 X 

5.1 X 

6.7 X 

5.8 X 



Table 6: Simulation results and theoretical values of nucleation rates for supersaturated methane, ethane, and carbon 
dioxide; bold values: threshold t > according to theory 



According to CNT, the size of the critical nucleus for 
ethane at 280 K, 2.80 mol/1, and 3397 kPa is = 78, 
cf. Table [6l The overheating of the critical nucleus as 
predicted by CNT from Equation (|8|) in this case, given 
that AHy = 8.3 kJ/mol and Z = 0.0089, amounts to 
AT* = 1.4 K; this corresponds to a nucleus temperature 
of 281.4 K that is actually observed for large nuclei in 
the simulation, but not for i ^ 78, cf. Figure [3] (center). 
From LFK we obtain ^* = 207 as weh as Z = 0.0062 and 
AT* = 0.98 K, which agrees well with the overheating ob- 
served for nuclei containing t ^ 207 molecules after a delay 
of 700 ps. 

Values of J{i) determined with the method of Yasuoka 
and Matsumoto are only significant for t > i* . Since the 
size of the critical nucleus can not be obtained by means 
of this method, the theories are checked against their own 
predictions of t*: if the theoretical value of l* is smaller 
than the threshold used to evaluate the MD simulation. 



then simulation and theory should be expected to agree. 
Such data are directly comparable - they correspond to 
the highlighted values in Tables [3l [H [5l and [6l 

The values collected for the quadrupolar fluids show 
an eccellent agreement with CNT: all directly compara- 
ble nucleation rates agree within one order of magnitude. 
For methane, CNT significantly underestimates the nucle- 
ation rate at 106 K and overestimates it at 170 K. The 
predictions of J based on the LFK model are generally 
too low for carbon dioxide, with an error between one and 
two orders of magnitude in the directly comparable cases. 
However, LFK predicts the nucleation rate accurately for 
methane at 106, 114, and 130 K as well as for ethane at 
176.5 K. 

Both theories are observed to deviate by about three 
orders of magnitude from certain directly comparable J{i) 
values: for methane at 106 K and 503 kPa, the method of 
Yasuoka and Matsumoto yields J(25) = 1.8 x 10^^/(m^s), 
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while CNT predicts J = 1.6 x 102^/(m^s), cf. TableU For 
methane at 170 K and Sp = 1.247, cf. Table [3l the nucle- 
ation rate according to the LFK model is J = 6.7 x 10^^/ 
(m^s) as opposed to J(225) = 2.7 x 10^^/(m^s), obtained 
from an MD simulation. 

Conclusion. Molecular dynamics simulations of large 
nucleating systems were conducted in the canonical en- 
semble and analyzed according to the method of Yasuoka 
and Matsumoto. It was shown that the average nu- 
cleus temperature, for a given nucleus size, decreases over 
time, which leads to a considerable increase of the net 
growth rates during simulation. Nucleation rate isotherms 



were obtained and compared to theoretical predictions for 
methane at 106, 114, 130, and 170 K, for ethane at 176.5 
and 280 K, as well as for carbon dioxide at 237, 253, 269, 
and 285 K. Nucleation rates from both CNT and LFK were 
observed to agree within three orders of magnitude with 
those simulation results that can directly be compared to 
the theories. In particular, CNT shows very small devia- 
tions for ethane and carbon dioxide over the entire temper- 
ature range. The LFK model consistently and significantly 
underpredicts nucleation rates at high temperatures, but 
agrees very well for methane and ethane at low tempera- 
tures. 




500 600 800 900 1000 1100 2000 2500 3500 4000 4500 500 1000 1500 

pressure / kPa pressure / kPa cluster size 



Figure 4: left - Nucleation rate of CIl4 at 106, 114, and 130 K; small squares: J(25), large squares: J(75) and J(150), 
circles: J(225) and J(700); center - Nucleation rate of CO2 at 237 and 269 K; small squares: J(25) and J(50), large 
squares: J(75), circles: J(300); right - First passage time of clusters in CH4 at 170 K with N = 250000 
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